In [1]:
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
from sklearn.utils import resample
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import scale
import plotly.io as pio
import plotly.graph_objects as go

import numpy as np
pd.options.display.float_format = '{:.5f}'.format
pio.renderers.default='notebook'
In [2]:
bd1 = pd.read_excel("./Datos/Base_datos_1.xlsx")
bd_V_farnesiana = bd1.loc[bd1['Especie'] == 'V. farnesiana'].drop(['Tejido', 'Replica', 'Especie','PESO '],axis=1)
bd_R_communis = bd1.loc[bd1['Especie'] == 'R. communis'].drop(['Tejido', 'Replica', 'Especie','PESO '],axis=1)
In [3]:
print('----------------------')
print('Media de cada variable')
print('----------------------')
bd1.groupby(["Especie","Trat"]).mean()
----------------------
Media de cada variable
----------------------
Out[3]:
Replica As Ca Cu Fe K Mg Mn P Pb Zn Biomasa
Especie Trat
R. communis C 12.00000 5.43667 16350.00000 84.03333 5585.00000 25616.66667 5210.00000 65.71667 2958.33333 6.48000 28.70000 5.48500
CI 73.50000 4.86625 16462.50000 94.15000 4720.00000 25600.00000 5785.00000 73.05000 3400.00000 6.68625 33.60000 5.96625
J 22.00000 10.32500 13900.00000 277.10000 14505.00000 13450.00000 5905.00000 58.70000 1335.00000 13.58000 40.70000 0.76500
JI 48.00000 15.60000 14550.00000 361.30000 18385.00000 10900.00000 4395.00000 64.40000 1335.00000 18.68000 48.00000 0.30500
SN 96.33333 4.41333 15683.33333 29.23000 1588.33333 11650.00000 5820.00000 84.03333 1915.00000 0.84500 88.03333 4.24333
V. farnesiana C 7.66667 6.31167 17850.00000 88.81167 4667.50000 14516.66667 2973.33333 93.60000 2210.00000 7.21000 36.06667 4.08500
CI 67.50000 8.01125 17025.00000 125.70500 5228.75000 18650.00000 3250.00000 95.01250 2095.00000 7.46000 31.61250 6.43000
J 28.00000 8.45000 12100.00000 245.55000 9825.00000 15550.00000 2550.00000 76.45000 965.00000 10.62500 51.20000 0.08545
JI 42.00000 9.93000 12500.00000 326.30000 9960.00000 14050.00000 2670.00000 108.80000 935.00000 9.84000 57.95000 0.23000
SN 89.00000 6.81750 11612.50000 57.21250 948.25000 13862.50000 2267.50000 84.01250 1831.25000 0.54000 69.06250 5.40125

Varianza de cada variable

La varianza es muy distinta entre las variables

In [4]:
bd1.groupby(["Especie","Trat"]).var()
Out[4]:
Replica As Ca Cu Fe K Mg Mn P Pb Zn Biomasa
Especie Trat
R. communis C 31.20000 23.27819 13619000.00000 5879.27367 25086230.00000 36585666.66667 11064720.00000 332.07767 174096.66667 23.01456 21.99200 2.35571
CI 28.85714 9.44388 9654107.14286 4835.82571 9040857.14286 30811428.57143 15135200.00000 242.19429 129657.14286 11.36663 34.59714 4.14043
J 0.00000 173.91125 8820000.00000 135096.02000 386142050.00000 6125000.00000 29722050.00000 480.50000 31250.00000 298.65680 134.48000 0.02645
JI 0.00000 414.72000 14045000.00000 232152.98000 620576450.00000 180000.00000 7411250.00000 2492.18000 110450.00000 586.18880 317.52000 0.00005
SN 8.26667 2.47227 11053666.66667 492.59404 1607176.66667 2723000.00000 14672120.00000 666.12267 127630.00000 0.17047 2827.86267 0.99683
V. farnesiana C 19.46667 43.78094 94547000.00000 8377.65554 24953595.10000 6485666.66667 3339386.66667 4432.36400 710160.00000 53.27256 796.57467 0.98035
CI 41.42857 69.51336 49193571.42857 16318.17054 32788207.64286 37094285.71429 4562200.00000 3008.32982 382885.71429 64.52351 471.68696 5.03086
J 0.00000 120.12500 62720000.00000 82783.80500 175781250.00000 26645000.00000 1548800.00000 400.44500 2450.00000 102.96125 74.42000 0.00130
JI 0.00000 152.07680 74420000.00000 152020.98000 154528200.00000 3645000.00000 217800.00000 2339.28000 2450.00000 133.17120 53.04500 0.01280
SN 56.57143 47.19822 4306964.28571 3066.85839 1032021.64286 5391250.00000 956135.71429 1517.76411 141955.35714 0.13954 4185.20839 3.09110

PCA Farnesiana

In [5]:
tratamientos = ["C","CI","J","JI","SN"]
for tratamiento in tratamientos:
    mi_db = bd_V_farnesiana.loc[bd_V_farnesiana['Trat'] == tratamiento].copy().drop(["Trat","Biomasa"],axis=1)
    mi_db2 = bd_V_farnesiana.loc[bd_V_farnesiana['Trat'] == tratamiento].copy().drop(["Trat"],axis=1)
    pca_pipe = make_pipeline(StandardScaler(), PCA())

    # Se entrena y extrae el modelo entrenado del pipeline
    modelo_pca = pca_pipe.fit(mi_db).named_steps['pca']

    # Se combierte el array a dataframe para aƱadir nombres a los ejes.
    componentes = []

    for componente in range(len(modelo_pca.components_)):
        nombre = "PC{}".format(componente+1)
        componentes.append(nombre)

    pca_df = pd.DataFrame(
        data    = modelo_pca.components_,
        columns = mi_db.columns,
        index   = componentes
    )

    #print(np.transpose(modelo_pca.components_))
    #print(modelo_pca.components_)
    pca_df_componentes = pd.DataFrame(
        data    = np.transpose(modelo_pca.components_),
        columns = componentes,
    )

    pca_df_variance = pd.DataFrame(
        data    = [modelo_pca.explained_variance_ratio_],
        columns = componentes,
    )

    prop_varianza_acum = modelo_pca.explained_variance_ratio_.cumsum()
    pca_df_acum = pd.DataFrame(
        data    = [prop_varianza_acum],
        columns = componentes,
    )
    df_acumpc = pca_df_acum.T
    df_acumpc.rename(columns = {0:"Value"}, inplace=True)

    #heatmap
    fig = px.imshow(pca_df)
    fig.update_layout(
        template="plotly_dark",
        title="Influencia de las variables en cada componente, Tratamiento: '%s'" % tratamiento,
    )
    fig.update_xaxes(title_text='Mineral')
    fig.update_yaxes(title_text='Componente')
    fig.show()

    #Variaza y varianza acumulada
    df_variance_pca = pca_df_variance.T
    df_variance_pca.rename(columns = {0:"Value"}, inplace=True)
    fig = go.Figure()
    fig.add_trace(
        go.Scatter(
            x=df_acumpc.index,
            y=df_acumpc['Value'],
            name = "Varianza Acumulada"
        ))
    fig.add_trace(
        go.Bar(
            x=df_variance_pca.index,
            y=df_variance_pca['Value'],
            name = "Varianza"
        ))
    fig.update_xaxes(title_text='Componente')
    fig.update_yaxes(
        range=(0, 1.1),
        constrain='domain',
        title_text='Varianza'
    )
    fig.update_layout(
        template="plotly_dark",
          title="Variaza y Varianza Acumulada, Tratamiento: '%s'" % tratamiento,
          legend=dict(
            orientation="h",
            yanchor="bottom",
            y=1.02,
            xanchor="right",
            x=1
        )
    )
    fig.show()

    proyecciones = pca_pipe.transform(X=mi_db)
    proyecciones = pd.DataFrame(
        proyecciones,
        columns = componentes,
        index   = mi_db.index
    )

    if(tratamiento=="J" or tratamiento=="JI"):
        fig = px.scatter(
            proyecciones, x="PC1", y="PC2",
            #title=f'Total Explained Variance: {total_var:.2f}%',
            labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'}
        )
        fig.update_layout(template="plotly_dark", title="Visualización de dos componentes, Tratamiento: '%s'" % tratamiento)
        fig.show()
    else:
        fig = px.scatter_3d(
            proyecciones, x="PC1", y="PC2", z="PC3", color=mi_db2["Biomasa"],
            #title=f'Total Explained Variance: {total_var:.2f}%',
            labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'}
        )
        fig.update_layout(
            template="plotly_dark",
            title="Visualización de tres componentes, Tratamiento: '%s'" % tratamiento, 
        )
        fig.show()

PCA Communis

In [6]:
tratamientos = ["C","CI","J","JI","SN"]
for tratamiento in tratamientos:
    mi_db = bd_R_communis.loc[bd_R_communis['Trat'] == tratamiento].copy().drop(["Trat","Biomasa"],axis=1)
    mi_db2 = bd_R_communis.loc[bd_R_communis['Trat'] == tratamiento].copy().drop(["Trat"],axis=1)
    pca_pipe = make_pipeline(StandardScaler(), PCA())

    # Se entrena y extrae el modelo entrenado del pipeline
    modelo_pca = pca_pipe.fit(mi_db).named_steps['pca']

    # Se combierte el array a dataframe para aƱadir nombres a los ejes.
    componentes = []

    for componente in range(len(modelo_pca.components_)):
        nombre = "PC{}".format(componente+1)
        componentes.append(nombre)

    pca_df = pd.DataFrame(
        data    = modelo_pca.components_,
        columns = mi_db.columns,
        index   = componentes
    )

    #print(np.transpose(modelo_pca.components_))
    #print(modelo_pca.components_)
    pca_df_componentes = pd.DataFrame(
        data    = np.transpose(modelo_pca.components_),
        columns = componentes,
    )

    pca_df_variance = pd.DataFrame(
        data    = [modelo_pca.explained_variance_ratio_],
        columns = componentes,
    )

    prop_varianza_acum = modelo_pca.explained_variance_ratio_.cumsum()
    pca_df_acum = pd.DataFrame(
        data    = [prop_varianza_acum],
        columns = componentes,
    )
    df_acumpc = pca_df_acum.T
    df_acumpc.rename(columns = {0:"Value"}, inplace=True)

    #heatmap
    fig = px.imshow(pca_df)
    fig.update_layout(
        template="plotly_dark",
        title="Influencia de las variables en cada componente, Tratamiento: '%s'" % tratamiento,
    )
    fig.update_xaxes(title_text='Mineral')
    fig.update_yaxes(title_text='Componente')
    fig.show()

    #Variaza y varianza acumulada
    df_variance_pca = pca_df_variance.T
    df_variance_pca.rename(columns = {0:"Value"}, inplace=True)
    fig = go.Figure()
    fig.add_trace(
        go.Scatter(
            x=df_acumpc.index,
            y=df_acumpc['Value'],
            name = "Varianza Acumulada"
        ))
    fig.add_trace(
        go.Bar(
            x=df_variance_pca.index,
            y=df_variance_pca['Value'],
            name = "Varianza"
        ))
    fig.update_xaxes(title_text='Componente')
    fig.update_yaxes(
        range=(0, 1.1),
        constrain='domain',
        title_text='Varianza'
    )
    fig.update_layout(
        template="plotly_dark",
          title="Variaza y Varianza Acumulada, Tratamiento: '%s'" % tratamiento,
          legend=dict(
            orientation="h",
            yanchor="bottom",
            y=1.02,
            xanchor="right",
            x=1
        )
    )
    fig.show()

    proyecciones = pca_pipe.transform(X=mi_db)
    proyecciones = pd.DataFrame(
        proyecciones,
        columns = componentes,
        index   = mi_db.index
    )

    if(tratamiento=="J" or tratamiento=="JI"):
        fig = px.scatter(
            proyecciones, x="PC1", y="PC2",
            #title=f'Total Explained Variance: {total_var:.2f}%',
            labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'}
        )
        fig.update_layout(template="plotly_dark", title="Visualización de dos componentes, Tratamiento: '%s'" % tratamiento)
        fig.show()
    else:
        fig = px.scatter_3d(
            proyecciones, x="PC1", y="PC2", z="PC3", color=mi_db2["Biomasa"],
            #title=f'Total Explained Variance: {total_var:.2f}%',
            labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'}
        )
        fig.update_layout(
            template="plotly_dark",
            title="Visualización de tres componentes, Tratamiento: '%s'" % tratamiento, 
        )
        fig.show()

PCA con Boostrap de los datos Farnesiana

In [7]:
n_iterations = 1000
new_data_df = pd.DataFrame()

tratamientos = ["C","CI","J","JI","SN"]

for tratamiento in tratamientos:
    mi_db = bd_V_farnesiana.loc[bd_V_farnesiana['Trat'] == tratamiento].copy().drop(["Trat"],axis=1)
    new_data_df = pd.DataFrame()
    #Boostrap
    for i in range(n_iterations):
        newData = resample(mi_db, replace=True, n_samples=len(mi_db))
        new_data_df = new_data_df.append(newData.mean().to_frame().T, ignore_index = True)
    new_data_df = new_data_df.append(mi_db, ignore_index = True)

    #PCA
    mi_db = new_data_df.drop(["Biomasa"],axis=1)
    mi_db2 = new_data_df.copy()
    pca_pipe = make_pipeline(StandardScaler(), PCA())

    # Se entrena y extrae el modelo entrenado del pipeline
    modelo_pca = pca_pipe.fit(mi_db).named_steps['pca']

    # Se combierte el array a dataframe para aƱadir nombres a los ejes.
    componentes = []

    for componente in range(len(modelo_pca.components_)):
        nombre = "PC{}".format(componente+1)
        componentes.append(nombre)

    pca_df = pd.DataFrame(
        data    = modelo_pca.components_,
        columns = mi_db.columns,
        index   = componentes
    )

    #print(np.transpose(modelo_pca.components_))
    #print(modelo_pca.components_)
    pca_df_componentes = pd.DataFrame(
        data    = np.transpose(modelo_pca.components_),
        columns = componentes,
    )

    pca_df_variance = pd.DataFrame(
        data    = [modelo_pca.explained_variance_ratio_],
        columns = componentes,
    )

    prop_varianza_acum = modelo_pca.explained_variance_ratio_.cumsum()
    pca_df_acum = pd.DataFrame(
        data    = [prop_varianza_acum],
        columns = componentes,
    )
    df_acumpc = pca_df_acum.T
    df_acumpc.rename(columns = {0:"Value"}, inplace=True)

    #heatmap
    fig = px.imshow(pca_df)
    fig.update_layout(
        template="plotly_dark",
        title="Influencia de las variables en cada componente, Tratamiento: '%s'" % tratamiento,
    )
    fig.update_xaxes(title_text='Mineral')
    fig.update_yaxes(title_text='Componente')
    fig.show()

    #Variaza y varianza acumulada
    df_variance_pca = pca_df_variance.T
    df_variance_pca.rename(columns = {0:"Value"}, inplace=True)
    fig = go.Figure()
    fig.add_trace(
        go.Scatter(
            x=df_acumpc.index,
            y=df_acumpc['Value'],
            name = "Varianza Acumulada"
        ))
    fig.add_trace(
        go.Bar(
            x=df_variance_pca.index,
            y=df_variance_pca['Value'],
            name = "Varianza"
        ))
    fig.update_xaxes(title_text='Componente')
    fig.update_yaxes(
        range=(0, 1.1),
        constrain='domain',
        title_text='Varianza'
    )
    fig.update_layout(
        template="plotly_dark",
          title="Variaza y Varianza Acumulada, Tratamiento: '%s'" % tratamiento,
          legend=dict(
            orientation="h",
            yanchor="bottom",
            y=1.02,
            xanchor="right",
            x=1
        )
    )
    fig.show()

    proyecciones = pca_pipe.transform(X=mi_db)
    proyecciones = pd.DataFrame(
        proyecciones,
        columns = componentes,
        index   = mi_db.index
    )

    if(tratamiento=="J" or tratamiento=="JI"):
        fig = px.scatter(
            proyecciones, x="PC1", y="PC2",
            #title=f'Total Explained Variance: {total_var:.2f}%',
            labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'}
        )
        fig.update_layout(template="plotly_dark", title="Visualización de dos componentes, Tratamiento: '%s'" % tratamiento)
        fig.show()
    else:
        fig = px.scatter_3d(
            proyecciones, x="PC1", y="PC2", z="PC3", color=mi_db2["Biomasa"],
            #title=f'Total Explained Variance: {total_var:.2f}%',
            labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'}
        )
        fig.update_layout(
            template="plotly_dark",
            title="Visualización de tres componentes, Tratamiento: '%s'" % tratamiento, 
        )
        fig.show()

PCA con Boostrap de los datos Communis

In [8]:
n_iterations = 1000
new_data_df = pd.DataFrame()

tratamientos = ["C","CI","J","JI","SN"]

for tratamiento in tratamientos:
    mi_db = bd_R_communis.loc[bd_R_communis['Trat'] == tratamiento].copy().drop(["Trat"],axis=1)
    new_data_df = pd.DataFrame()
    #Boostrap
    for i in range(n_iterations):
        newData = resample(mi_db, replace=True, n_samples=len(mi_db))
        new_data_df = new_data_df.append(newData.mean().to_frame().T, ignore_index = True)
    new_data_df = new_data_df.append(mi_db, ignore_index = True)

    #PCA
    mi_db = new_data_df.drop(["Biomasa"],axis=1)
    mi_db2 = new_data_df.copy()
    pca_pipe = make_pipeline(StandardScaler(), PCA())

    # Se entrena y extrae el modelo entrenado del pipeline
    modelo_pca = pca_pipe.fit(mi_db).named_steps['pca']

    # Se combierte el array a dataframe para aƱadir nombres a los ejes.
    componentes = []

    for componente in range(len(modelo_pca.components_)):
        nombre = "PC{}".format(componente+1)
        componentes.append(nombre)

    pca_df = pd.DataFrame(
        data    = modelo_pca.components_,
        columns = mi_db.columns,
        index   = componentes
    )

    #print(np.transpose(modelo_pca.components_))
    #print(modelo_pca.components_)
    pca_df_componentes = pd.DataFrame(
        data    = np.transpose(modelo_pca.components_),
        columns = componentes,
    )

    pca_df_variance = pd.DataFrame(
        data    = [modelo_pca.explained_variance_ratio_],
        columns = componentes,
    )

    prop_varianza_acum = modelo_pca.explained_variance_ratio_.cumsum()
    pca_df_acum = pd.DataFrame(
        data    = [prop_varianza_acum],
        columns = componentes,
    )
    df_acumpc = pca_df_acum.T
    df_acumpc.rename(columns = {0:"Value"}, inplace=True)

    #heatmap
    fig = px.imshow(pca_df)
    fig.update_layout(
        template="plotly_dark",
        title="Influencia de las variables en cada componente, Tratamiento: '%s'" % tratamiento,
    )
    fig.update_xaxes(title_text='Mineral')
    fig.update_yaxes(title_text='Componente')
    fig.show()

    #Variaza y varianza acumulada
    df_variance_pca = pca_df_variance.T
    df_variance_pca.rename(columns = {0:"Value"}, inplace=True)
    fig = go.Figure()
    fig.add_trace(
        go.Scatter(
            x=df_acumpc.index,
            y=df_acumpc['Value'],
            name = "Varianza Acumulada"
        ))
    fig.add_trace(
        go.Bar(
            x=df_variance_pca.index,
            y=df_variance_pca['Value'],
            name = "Varianza"
        ))
    fig.update_xaxes(title_text='Componente')
    fig.update_yaxes(
        range=(0, 1.1),
        constrain='domain',
        title_text='Varianza'
    )
    fig.update_layout(
        template="plotly_dark",
          title="Variaza y Varianza Acumulada, Tratamiento: '%s'" % tratamiento,
          legend=dict(
            orientation="h",
            yanchor="bottom",
            y=1.02,
            xanchor="right",
            x=1
        )
    )
    fig.show()

    proyecciones = pca_pipe.transform(X=mi_db)
    proyecciones = pd.DataFrame(
        proyecciones,
        columns = componentes,
        index   = mi_db.index
    )

    if(tratamiento=="J" or tratamiento=="JI"):
        fig = px.scatter(
            proyecciones, x="PC1", y="PC2",
            #title=f'Total Explained Variance: {total_var:.2f}%',
            labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'}
        )
        fig.update_layout(template="plotly_dark", title="Visualización de dos componentes, Tratamiento: '%s'" % tratamiento)
        fig.show()
    else:
        fig = px.scatter_3d(
            proyecciones, x="PC1", y="PC2", z="PC3", color=mi_db2["Biomasa"],
            #title=f'Total Explained Variance: {total_var:.2f}%',
            labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'}
        )
        fig.update_layout(
            template="plotly_dark",
            title="Visualización de tres componentes, Tratamiento: '%s'" % tratamiento, 
        )
        fig.show()